Pre-visualization of the normalized count data before any differential analysis

Table of content

  1. Sample clustering
    1. With all genes
    2. Without X/Y genes
  2. Check samples and genes with normalized counts
    1. Check for artefact due to FACS sorting
    2. Check for contamination from other cell types that could have escaped the sorting-gating
  3. PCA on the normalized counts
    1. With all factors
    2. SPF vs GF
  4. Gene co-expression analysis
    1. Extract modules of co-expressed genes
    2. Relationship between modules and samples
    3. Genes in modules
      1. Heatmaps with Z-scores
      2. Sinaplots of the Z-scores per groups
      3. Median of counts over ages
    4. Enrichment analysis in modules
    5. Network visualization
    6. Functional gene network

Load

Libraries and functions

Warning message in is.na(x[[i]]):
“is.na() applied to non-(list or vector) of type 'environment'”Warning message in rsqlite_fetch(res@ptr, n = n):
“Don't need to call dbFetch() for statements, only for queries”
==========================================================================
*
*  Package WGCNA 1.63 loaded.
*
*    Important note: It appears that your system supports multi-threading,
*    but it is not enabled within WGCNA in R. 
*    To allow multi-threading within WGCNA with all available cores, use 
*
*          allowWGCNAThreads()
*
*    within R. Use disableWGCNAThreads() to disable threading if necessary.
*    Alternatively, set the following environment variable on your system:
*
*          ALLOW_WGCNA_THREADS=<number_of_processors>
*
*    for example 
*
*          ALLOW_WGCNA_THREADS=4
*
*    To set the environment variable in linux bash shell, type 
*
*           export ALLOW_WGCNA_THREADS=4
*
*     before running R. Other operating systems or shells will
*     have a similar command to achieve the same aim.
*
==========================================================================


Allowing multi-threading with up to 4 threads.
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."

Data

Check samples and specific genes

Check for artefact due to FACS sorting

List of genes to check: Egr1, Jun, Zfp36l1, Malat1, Dusp1, Nr4a1, Fos

Normalized counts

Z-scores

Check for contamination from other cell types that could have escaped the sorting gating

Warning message:
“Setting row names on a tibble is deprecated.”Warning message:
“Setting row names on a tibble is deprecated.”Warning message:
“Setting row names on a tibble is deprecated.”

Remove genes not in list

cell
Scl2a5microglia
Ly6ggranulocyte
Trbc1t
Trbc2t
Igkcb
Ighmb
Mcpt4mast
Cma1mast
NKp46nk
NK1.1nk
NKG2Dnk
Tbetnk

Normalized counts

Z-score

Filter contamination genes and contaminated samples

Genes that could be artefact due to FACS sorting or contamination from other cell types and samples that were potentially contaminated

  1. 19694
  2. 61
  1. 19662
  2. 58
  1. 58
  2. 4

Sample clustering

With all genes

Without X/Y genes

PCA on the normalized counts

With all factors

SPF vs GF

M, 8w

M, 52w

M, 104w

F, 8w

F, 52w

F, 104w

Gene co-expression analysis

Extract modules of co-expressed genes

Filter genes that have a count > 10 in more than 90% of the samples:

  1. 19694
  2. 61
  1. 10291
  2. 61
  1. 10277
  2. 58
pickSoftThreshold: will use block size 4353.
 pickSoftThreshold: calculating connectivity for given powers...
   ..working on genes 1 through 4353 of 10277
   ..working on genes 4354 through 8706 of 10277
   ..working on genes 8707 through 10277 of 10277
   Power SFT.R.sq  slope truncated.R.sq  mean.k. median.k. max.k.
1      1    0.116  0.852          0.943 1040.000  1.04e+03 1630.0
2      2    0.353 -1.100          0.941  336.000  3.18e+02  821.0
3      3    0.664 -1.650          0.978  136.000  1.17e+02  496.0
4      4    0.769 -1.920          0.983   64.300  4.86e+01  332.0
5      5    0.806 -2.060          0.979   33.700  2.18e+01  238.0
6      6    0.813 -2.140          0.969   19.100  1.04e+01  178.0
7      7    0.817 -2.190          0.965   11.600  5.21e+00  137.0
8      8    0.842 -2.150          0.975    7.330  2.71e+00  108.0
9      9    0.864 -2.100          0.983    4.840  1.47e+00   87.5
10    10    0.876 -2.060          0.986    3.310  8.04e-01   71.6
11    12    0.897 -1.990          0.991    1.670  2.64e-01   49.7
12    14    0.904 -1.930          0.991    0.921  9.27e-02   35.7
13    16    0.919 -1.870          0.994    0.541  3.43e-02   26.4
14    18    0.930 -1.810          0.997    0.335  1.34e-02   19.9
15    20    0.940 -1.760          0.998    0.216  5.39e-03   15.2
PowerSFT.R.sqslopetruncated.R.sqmean.k.median.k.max.k.
1 0.1160946 0.8516606 0.9432081 1043.23197431.039858e+031631.43013
2 0.3534325 -1.0957144 0.9406479 335.81363833.181957e+02 820.65787
3 0.6637248 -1.6460594 0.9780977 136.38819671.174060e+02 496.21168
4 0.7685466 -1.9220353 0.9832291 64.31907084.862867e+01 332.43621
5 0.8057060 -2.0564238 0.9792358 33.71886882.180012e+01 237.66372
6 0.8134695 -2.1404068 0.9685734 19.13648601.041504e+01 177.64476
7 0.8167563 -2.1947647 0.9653164 11.55236995.214321e+00 137.13604
8 0.8422749 -2.1489288 0.9752509 7.32730562.709162e+00 108.47588
9 0.8638155 -2.1041594 0.9828091 4.83919561.465842e+00 87.45249
10 0.8760665 -2.0646700 0.9856232 3.30534528.042653e-01 71.58565
12 0.8974660 -1.9875884 0.9909974 1.67249372.641405e-01 49.69122
14 0.9036729 -1.9326232 0.9905956 0.92116289.265626e-02 35.74188
16 0.9189432 -1.8686196 0.9937579 0.54149043.432230e-02 26.39679
18 0.9297761 -1.8072410 0.9965048 0.33515971.342172e-02 19.89803
20 0.9401800 -1.7582833 0.9982447 0.21632965.386947e-03 15.24642
Powermean.k.
5 5 33.7188688
6 6 19.1364860
7 7 11.5523699
8 8 7.3273056
9 9 4.8391956
1010 3.3053452
1112 1.6724937
1214 0.9211628
1316 0.5414904
1418 0.3351597
1520 0.2163296

Parameters for network construction

 Calculating module eigengenes block-wise from all genes
   Flagging genes and samples with too many missing values...
    ..step 1
 ..Working on block 1 .
    TOM calculation: adjacency..
    ..will use 4 parallel threads.
     Fraction of slow calculations: 0.000000
    ..connectivity..
    ..matrix multiplication (system BLAS)..
    ..normalization..
    ..done.
   ..saving TOM for block 1 into file norm_genes_TOM-block.1.RData
 ....clustering..
 ....detecting modules..
 ....calculating module eigengenes..
 ....checking kME in modules..
     ..removing 45 genes from module 1 because their KME is too low.
     ..removing 27 genes from module 2 because their KME is too low.
     ..removing 19 genes from module 3 because their KME is too low.
     ..removing 7 genes from module 4 because their KME is too low.
     ..removing 1 genes from module 5 because their KME is too low.
     ..removing 2 genes from module 6 because their KME is too low.
     ..removing 4 genes from module 7 because their KME is too low.
     ..removing 4 genes from module 8 because their KME is too low.
     ..removing 1 genes from module 9 because their KME is too low.
     ..removing 1 genes from module 10 because their KME is too low.
     ..removing 1 genes from module 11 because their KME is too low.
     ..removing 1 genes from module 14 because their KME is too low.
     ..removing 1 genes from module 15 because their KME is too low.
     ..removing 1 genes from module 16 because their KME is too low.
     ..removing 2 genes from module 18 because their KME is too low.
     ..removing 1 genes from module 24 because their KME is too low.
 ..merging modules that are too close..
     mergeCloseModules: Merging modules whose distance is less than 0.35
       Calculating new MEs...

 ME0  ME1  ME2  ME3  ME4  ME5  ME6  ME7  ME8  ME9 
2719 1627 1240 1221 1123  925  467  467  371  117
9
  • 9 modules with between 117 and 1627 genes
  • 2719 genes not assigned to a module
7558
0
'grey'
1
'turquoise'
2
'blue'
3
'brown'
4
'yellow'
5
'green'
6
'red'
7
'black'
8
'pink'
9
'magenta'

Relationship between modules and samples

pdf: 2

Genes in modules

Percentages of XY genes per modules

[1] "0: 83 - 3.05259286502391%"
[1] "1: 59 - 3.62630608481868%"
[1] "2: 30 - 2.41935483870968%"
[1] "3: 58 - 4.75020475020475%"
[1] "4: 26 - 2.31522707034728%"
[1] "5: 27 - 2.91891891891892%"
[1] "6: 7 - 1.49892933618844%"
[1] "7: 12 - 2.56959314775161%"
[1] "8: 7 - 1.88679245283019%"
[1] "9: 6 - 5.12820512820513%"

Heatmaps with Z-scores

Type / gender / age

Gender / type / age

Age / type / gender

Sinaplots of the Z-scores per groups

The mean of the Z-score over the samples in the group is plot for each gene

Type / Age

Type / Age / Gender

Age / Type / Gender

Median of counts over ages

Enrichment analysis in modules

[1] "0"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
11858 GO:0044237            1.233688e-18                        1       1477
11875 GO:0044260            2.358207e-18                        1       1243
3933  GO:0008152            2.445785e-17                        1       1588
9635  GO:0034660            6.185192e-17                        1        124
11359 GO:0043170            2.606555e-16                        1       1320
16801 GO:0071704            1.109589e-15                        1       1522
      numInCat                                     term ontology
11858     8369               cellular metabolic process       BP
11875     6828 cellular macromolecule metabolic process       BP
3933      9196                        metabolic process       BP
9635       405                  ncRNA metabolic process       BP
11359     7434          macromolecule metabolic process       BP
16801     8810      organic substance metabolic process       BP
[1] "1"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
17654 GO:0090304            7.170038e-22                        1        568
11875 GO:0044260            2.938406e-19                        1        818
5709  GO:0016070            3.571231e-19                        1        512
2982  GO:0006139            2.636737e-17                        1        590
12846 GO:0046483            3.842591e-16                        1        592
11359 GO:0043170            4.718493e-16                        1        852
      numInCat                                             term ontology
17654     4158                   nucleic acid metabolic process       BP
11875     6828         cellular macromolecule metabolic process       BP
5709      3717                            RNA metabolic process       BP
2982      4658 nucleobase-containing compound metabolic process       BP
12846     4751                    heterocycle metabolic process       BP
11359     7434                  macromolecule metabolic process       BP
[1] "2"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
11858 GO:0044237            2.323529e-31                        1        748
3933  GO:0008152            1.098594e-29                        1        800
16801 GO:0071704            1.487654e-29                        1        772
11859 GO:0044238            2.230477e-29                        1        739
3427  GO:0006807            2.610162e-22                        1        523
9621  GO:0034641            5.680856e-22                        1        498
      numInCat                                         term ontology
11858     8369                   cellular metabolic process       BP
3933      9196                            metabolic process       BP
16801     8810          organic substance metabolic process       BP
11859     8342                    primary metabolic process       BP
3427      5576          nitrogen compound metabolic process       BP
9621      5271 cellular nitrogen compound metabolic process       BP
[1] "3"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
5311  GO:0015031            7.840882e-14                        1        189
12253 GO:0045184            1.767297e-13                        1        202
11873 GO:0044257            3.526414e-12                        1         90
11878 GO:0044265            3.966909e-12                        1        106
11644 GO:0043632            1.355917e-11                        1         78
14337 GO:0051603            1.837011e-11                        1         85
      numInCat                                                       term
5311      1603                                          protein transport
12253     1762                      establishment of protein localization
11873      604                         cellular protein catabolic process
11878      769                   cellular macromolecule catabolic process
11644      497     modification-dependent macromolecule catabolic process
14337      572 proteolysis involved in cellular protein catabolic process
      ontology
5311        BP
12253       BP
11873       BP
11878       BP
11644       BP
14337       BP
[1] "4"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
17654 GO:0090304            1.043412e-24                        1        410
2982  GO:0006139            1.107236e-22                        1        436
12846 GO:0046483            2.450298e-22                        1        441
3382  GO:0006725            3.346118e-21                        1        441
9621  GO:0034641            2.281973e-20                        1        465
5709  GO:0016070            2.728075e-20                        1        363
      numInCat                                             term ontology
17654     4158                   nucleic acid metabolic process       BP
2982      4658 nucleobase-containing compound metabolic process       BP
12846     4751                    heterocycle metabolic process       BP
3382      4812     cellular aromatic compound metabolic process       BP
9621      5271     cellular nitrogen compound metabolic process       BP
5709      3717                            RNA metabolic process       BP
[1] "5"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
3933  GO:0008152            3.891028e-27                        1        578
11858 GO:0044237            5.508154e-27                        1        533
9621  GO:0034641            4.944689e-24                        1        374
3552  GO:0007005            3.487645e-22                        1         98
3427  GO:0006807            4.233611e-22                        1        384
11859 GO:0044238            5.250566e-21                        1        509
      numInCat                                         term ontology
3933      9196                            metabolic process       BP
11858     8369                   cellular metabolic process       BP
9621      5271 cellular nitrogen compound metabolic process       BP
3552       610                   mitochondrion organization       BP
3427      5576          nitrogen compound metabolic process       BP
11859     8342                    primary metabolic process       BP
[1] "6"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
5709  GO:0016070            3.486898e-12                        1        171
4789  GO:0010467            1.850682e-10                        1        184
17654 GO:0090304            4.849420e-10                        1        177
13827 GO:0050794            1.233021e-09                        1        306
13823 GO:0050789            1.450196e-09                        1        318
3096  GO:0006351            1.831491e-09                        1        140
      numInCat                             term ontology
5709      3717            RNA metabolic process       BP
4789      4448                  gene expression       BP
17654     4158   nucleic acid metabolic process       BP
13827     8769   regulation of cellular process       BP
13823     9254 regulation of biological process       BP
3096      2987     transcription, DNA-templated       BP
[1] "7"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
12052 GO:0044710            2.245565e-11                        1        144
2881  GO:0005975            1.628591e-08                        1         37
3430  GO:0006810            1.859964e-08                        1        154
18602 GO:1901135            2.504175e-08                        1         53
14117 GO:0051234            7.722787e-08                        1        156
16800 GO:0071702            1.890205e-07                        1         99
      numInCat                                      term ontology
12052     3276         single-organism metabolic process       BP
2881       489            carbohydrate metabolic process       BP
3430      3906                                 transport       BP
18602      911 carbohydrate derivative metabolic process       BP
14117     4053             establishment of localization       BP
16800     2269               organic substance transport       BP
[1] "8"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
3522  GO:0006955            1.811570e-25                        1         76
1045  GO:0002376            1.523749e-22                        1        106
3519  GO:0006952            8.144335e-19                        1         72
12200 GO:0045087            5.453585e-17                        1         44
11388 GO:0043207            2.645861e-15                        1         48
14385 GO:0051707            2.645861e-15                        1         48
      numInCat                                 term ontology
3522       974                      immune response       BP
1045      1958                immune system process       BP
3519      1148                     defense response       BP
12200      499               innate immune response       BP
11388      679 response to external biotic stimulus       BP
14385      679           response to other organism       BP
[1] "9"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
3139  GO:0006412            4.272001e-32                        1         44
11288 GO:0043043            1.373090e-31                        1         44
11625 GO:0043604            2.216043e-30                        1         44
3213  GO:0006518            6.237136e-29                        1         44
11624 GO:0043603            2.528913e-28                        1         45
18738 GO:1901566            1.500784e-23                        1         45
      numInCat                                         term ontology
3139       556                                  translation       BP
11288      575                 peptide biosynthetic process       BP
11625      639                   amide biosynthetic process       BP
3213       691                    peptide metabolic process       BP
11624      798             cellular amide metabolic process       BP
18738     1097 organonitrogen compound biosynthetic process       BP

Network visualization

TOM calculation: adjacency..
..will use 4 parallel threads.
 Fraction of slow calculations: 0.000000
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.
0610007P14Rik0610009B22Rik0610009O20Rik0610010K14RikMettl260610012G03Rik0610030E20Rik0610037L13RikFmc11110002L01Rik⋯Zw10ZwintZxdaZxdbZxdcZyg11bZyxZzef1Zzz3Hikeshi
0610007P14Rik1111001010⋯0100000111
0610009B22Rik1110000010⋯0111001001
0610009O20Rik1110000100⋯1111011011
0610010K14Rik1001111010⋯0010011110
Mettl260001110010⋯1010010000
0610012G03Rik0001110010⋯1010000000